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Computational aeroacoustics requires efficient, high-resolution simulation tools. For smooth problems, this is 
best accomplished with very high-order in space and time methods on small stencils. However, the complexity of 
highly accurate numerical methods can inhibit their practical application, especially in irregular geometries. This 
complexity is reduced by using a special form of Hermite divided-difference spatial interpolation on Cartesian 
grids, and a Cauchy-Kowalewski recursion procedure for time advancement. In addition, a stencil constraint tree 
reduces the complexity of interpolating grid points that are located near wall boundaries. These procedures are 
used to develop automatically and to implement very high-order methods (>15) for solving the linearized Euler 
equations that can achieve less than one grid point per wavelength resolution away from boundaries by including 
spatial derivatives of the primitive variables at each grid point The accuracy of stable surface treatments is 
currently limited to 11th order for grid aligned boundaries and to 2nd order for irregular boundaries. 


I. Introduction 

N OISE generation and propagation are difficult to simulate nu- 
merically for a variety of well -documented reasons and require 
high-order numerical schemes. 12 However, high-order schemes can 
introduce a number of complications, such as the following: 

1) Large stencils near boundaries, with either many fictitious grid 
points, or large one-sided stencils, introduce programming com- 
plexity and numerical instability. 3 4 

2) High-order finite difference equations can require boundary 
treatments beyond the physical conditions of the original problem, 
which can excite spurious waves and instabilities. 5 

3) Generation of high-order, smooth, body-fitted grids around 
complex configurations can be difficult. 6 

4) High-order formulations can lack nonlinear robustness. 6 

5) The general usefulness of high-order methods is limited by 
first-order accurate shock capturing. 7 

The approach described in this paper will focus on the first three 
issues. The fourth issue is being investigated. 8 The last issue may 
possibly be avoided by including the physical viscosity and resolv- 
ing the steep gradients, 9 but with very high-resolution methods. 

In addressing the first three concerns, we limit ourselves to wave 
propagation and scattering problems governed by the linearized Eu- 
ler equations. Previous studies point out the advantages of high- 
order methods for acoustical propagation. 10 ' 15 Many methods in 
general use stop at fourth-order accuracy for time-dependent prob- 
lems because they use Runge-Kutta methods. High-order Runge- 
Kutta methods become notoriously difficult to derive because the 
number of nonlinear order conditions that need to be solved grows 
exponentially, that is, a 1 2th-order method has 78 1 3 nonlinear order 
conditions. 16 21 The advantages of using Runge-Kutta methods at 
orders less than 6 are commonly cited as flexibility, large stability 
limits, and ease of programming. 22 The practical limit on their order 
has been an impediment to the analysis of their use in high-order 
approaches for time -dependent acoustic applications. 23 

In this paper we use a series of explicit, local, high-order methods 
that have the same order of accuracy in space as in time. 24,25 These 
methods use Hermite interpolation on stencils that are two points 
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wide and a Cauchy-Kowalewski recursive procedure 26 for obtaining 
time derivatives from the space derivatives of the interpolant. The 
Lime derivatives are then used to advance the primitive variables and 
their spatial derivatives in time with a Taylor series expansion. This 
general approach is called the modified expansion solution approx- 
imation (MESA) method. 11 This method can be used to derive and 
implement algorithms with arbitrarily high orders of accuracy in 
multiple space dimensions if their complexity is properly managed 
and the computer’s floating point precision is sufficiently high. 27 

The complex task of developing and coding a multidimensional 
interpolant for each MESA scheme can be eliminated using the 
tensor product 28 of a new divided-difference form of Hermitian 
spatial interpolation on a two-point stencil. 27 The task of obtain- 
ing time derivatives with Cauchy-Kowalewski recursion 26 can be 
implemented by unrolling the recursion. 29 With these techniques, 
very high-order MESA schemes may be implemented using only a 
few pages of code, and their accuracy may be adjusted by merely 
changing loop indices, not the code. 27 This in turn enables local 
solution adaptive order changes at each time step. 

Most unsteady flow simulators using Cartesian grids in complex 
geometries are based on a finite volume approach 30 with the excep- 
tion of the work by Kurbatskii and Tam, 31 which specifically in- 
vestigates acoustic scattering using a fourth-order accurate in time 
finite difference approach. In their work, curved surfaces are approx- 
imated with linear segments and ghost points are used to enforce 
boundary conditions, all of which incurred a fair amount of pro- 
gramming complexity to implement. 31 

The approach presented here, however, does not need to approxi- 
mate the geometry and the programming task is completed automat- 
ically using computer algebra 29 in a manner similar to the works of 
Wirth 32 and Steinberg and Roache. 33 This is accomplished by view- 
ing the boundary conditions as applying uniformly in time through- 
out the boundary surface. The boundary conditions can be differen- 
tiated in the surface and in time, and the governing equations can 
be used to obtain an infinite number of constraints in the boundary. 4 
A stencil constraint tree is used to simplify the task of symbolically 
imposing these high-order surface boundary conditions. 29 

With these techniques, the accuracy of interior propagation seems 
limited only by the floating point precision available. However, we 
stably attained only llth-order accuracy for acoustic scattering in 
a box with sides that are aligned to the grid, and we attained only 
2nd-order accuracy in more general cases. These limits are due to 
poorly conditioned matrix systems and numerical instability. How- 
ever, these effects are dependent on the choice of boundary condi- 
tions and require further study as will be shown. 

The objective of this paper is to present and validate a new ap- 
proach to aeroacoustic computations in complex geometries that 
has the potential to fully utilize the precision of today’s computers. 
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Fig. 1 Cascade superimposed on grid. 
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A common theme in this paper will be the many advantages of 
using a two-point wide stencil that include very high resolution, so- 
lution adaptability, and Courant-Friedrichs-Lewy (CFL) stability 
near boundaries. One of the more intriguing results to be presented 
here is the possibility of subgrid-scale resolution using solution 
adaptable algorithms. 

II. Governing Equations 

We will demonstrate this new approach by solving the linearized 
Euler equations in a uniform mean flowfield with mean velocity M , 
perturbation velocity u, pressure p, and Cartesian coordinates x t . 
We assume the initial conditions are known and that no additional 
sources are present. The conservation of momentum and energy 
equations are 


diij 3 Uj dp 

— - + M, — - H — = 0, 

3/ 3 Xt 3 xj 


dp dp 3 Uj ^ ... 

s 7 + m '^ + S7‘ 0 (1) 


and are nondimensionalized with respect to the following scales: 
length scale L, velocity scale c (sound speed), timescale L/c , density 
scale pQ (ambient density), and pressure scale p^c 2 . Also note that 
the continuity equation is not included because it is not necessary 
for calculating the acoustic response. 

We will use 2s + 1 -order explicit MESA methods, where s can be 
any nonnegative integer. We will use only two-point stencils, where 
each grid point contains the primitive solution variables p and w, 
and their spatial derivatives (Hermitian data), for a total of 3 (j + 1 ) 2 
or 4 (s + l) 3 terms per grid point in two or three spatial dimensions, 
respectively. For example, in two spatial dimensions with s = 1, we 
have a third-order method that uses the following data on each grid 
point: p , m, u, p x , u x , v x , p y , u yy v yy p xy , u xyy and v xy . 


III. Grid Classification 

Generally, we are interested in simulating the acoustical scattering 
from objects that are defined parametrically and that are superim- 
posed on a uniform Cartesian grid, as shown in Fig. 1 for a cascade 
of three airfoils. Notice that each grid point is labeled with a closed 
circle, an open circle, or the letter B, and represents an interior, a 
fill, or a boundary grid point, respectively. In Fig. 1, the assumed 
stencil size is three points per dimension, and a fill grid point is 
simply one in which one of its neighboring grid points is either on 
or beyond the boundary. It is referred to as a fill because that grid 
point cannot be computed directly and needs to be filled with data. 
If none of the neighboring grid points is on or beyond the boundary, 
then it is classified as an interior grid point and can be computed 
normally. All other grid points are considered boundary points and 
are not needed because ghost grid points are not used. 

Grid point classification proceeds by first labeling all grid points 
as boundary type and then finding the interior and fill -type grid 
points with a simple recursive procedure. The procedure determines 
if any neighboring grid points are boundary points by checking if 
the surface of any object intersects the imaginary line joining the 
center point to its neighbor. Because we are using a finite difference 
method, we need to check only for line-to-surface intersections, 
which is much simpler than determining the surface-to-surface in- 
tersections necessary if a finite volume method is used. If there is no 
intersection, then the current grid point is a fill point and recursion 
stops. Otherwise it is an interior point, and the procedure is called 
again, but starting with the neighboring grid points. 

IV. Advancing the Solution Using MESA 

With the uniform Cartesian grid points correctly labeled and the 
initial data assigned to all of the interior and fill grid points, the next 
step is to advance the solution in time. Because interior grid points 
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Fig. 2 Staggered grid using two-point stencil. 
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and their stencils are completely contained within the computational 
domain, they will be advanced directly using an efficient form of 
the MESA schemes. However, because a fill point’s stencil reaches 
outside the computational domain field, its data will be obtained 
by spatial interpolation. At each time step, every grid point is first 
treated as an interior point to simplify the looping structure for 
efficient vector and parallel operations, and then the relatively few 
fill points are recomputed using interpolation. 

Advancing the interior points with MESA is a three-step process. 
The spatial derivatives at the center of the stencil are approximated 
using spatial interpolation, the spatial derivatives are converted to 
time derivatives, and the solution is advanced in time with a Taylor 
series. The MESA scheme 34 is a general approach to developing 
high-order numerical schemes, but for this work we are only inter- 
ested in methods with stencils that are two grid points wide in all 
space dimensions because they enable arbitrarily high-order spa- 
tial interpolation with no decrease in CFL stability bounds near a 
boundary. 

These two-point stencil schemes are numerically stable only when 
a staggered grid is used at each time step 34 as shown in Fig. 2 for a 
two-dimensional stencil. The locations marked with Xs in Fig. 2 are 
advanced first using the 2 x 2 stencils centered about each X. Next, 
the Hermitian data now stored at the X locations are used as initial 
data to advance the solution to the center as indicated by the large 
dot in Fig. 2. Staggering the grid has the same effect as applying 
the MESA scheme to the entire 3x3 stencil in the Fig. 2 and then 
adding artificial dissipation. It is, however, more efficient to stagger 
the grid because neighboring grid points will reuse the data at the X 
locations in Fig. 2. In two dimensions, this staggered grid procedure 
is numerically stable if 34 

At 1 

X = — < (2) 

Ax ~ 1 + maxflA/jJ, iM v |} 


A. Hermitian Derivative Approximation 

The MESA scheme requires an approximation to the solution of 
the primitive variables p and u t and their spatial derivatives at the 
center of each stencil. Let the function f(x,y) be an approximation 
of one of the primitive variables from Eqs. (1) and define its origin 
to be at the center of each two-point stencil shown in Fig. 3. Then in 
two dimensions, the following data must be approximated once for 
each primitive variable [/(x, y) = p(x , y), u(x, y), and v(jt, y)]: 

S a ^ b f( 0, 0) 

J ' V a,b:a,b = 0,1,2 25+1 (3) 

dx a dy b 


On a Cartesian grid, these approximations can be found using 
tensor product interpolation 28 in two or three dimensions if the 
interpolant 


/(^)=EE c i^ 


with its origin at stencil center is used with a 25 + 1 -order method 
and if the following data is available at each grid point: 

d a ^ b f(x,y) 

3x a dy b w 

Determining the C[j = f '/(x, y)/3x i y j ] terms in 

the two-dimensional interpolant [Eq. (4)] for very high-order meth- 
ods is very difficult and inefficient unless it is reduced to a series 
of one-dimensional interpolations by using tensor products. Briefly, 
this is accomplished by performing one-dimensional interpolations 
of the form 27 


m = J2 c - x ' 


Fig. 3 Two-point two-dimensional 
stencil. 
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It is then possible to interpolate in the x direction, using the data 
on the grid points, the following two sets of data for y = y 0 and yj : 


(3- 


f(x=0, y) 
dx'dyJ 


i =0, 1,2, ..., 25 + 1 

V L J • . ni (7) 

J = 0, 1 5 


Then these data are used to interpolate in the y direction using 
the following one-dimensional shape function: 

Is - ] 

f(x=0,y) = J2 C fy J <*) 

> = 0 

Once the Cj terms are determined, they provide the data shown in 
Eq. (3): 

, / 1 \ 3- =>7(jc =0,y = 0) 

1 \i\j\) 3x'dyJ 

The one-dimensional shape functions [Eqs. (6) and (8)] can be di- 
rectly solved using computer algebra (Groebner bases) (see Ref. 1 8) 
to create algebraic expressions for each term for up to ap- 
proximately 30th-order methods. However, this results in lengthy, 
inefficient expressions that limit the accuracy and prevent instan- 
taneous accuracy changes necessary for resolving varying wave- 
length scales. A better way is to use Hermitian divided-differences, 
which will create the equivalent one -dimensional shape function 
(Newton’s interpolatory formula), 35 

2(f hJ-1) - 1 

fix) = 22 Qu(x-x 0 y + l ix-x + 


7^ Qi.iix - x 0 y 


where the Q iti are the traditional Hermitian divided differences from 
the tableau (shown in Fig. 4) which are given by 

DO i = 0, 5 


DO j = 0, i 




END DO 


END DO 


DO i = 5 + 1, 2*5+1 
DO j = i — 5, i 

„ Qi.i - 1 - Qi-i.j- 


ENDDO 


(6) 


END DO 
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Fig. 4 

Fifth-order divided -difference tableau. 


It is not enough to simply find Newton’s interpolator formula 
[Eq. (10)] because we need to approximate p , u , and v and their 
spatial derivatives at the stencil center. This would ordinarily require 
differentiating Newton’s Eq. (10) using computer algebra. By the 
product rule of differentiation, Newton’s formula will double in size 
after each differentiation, and a 50th-order MESA scheme would 
need a Newton’s formula of approximately 2 50 ~ 10 16 terms. This 
issue is eliminated by the following result, which is applicable only 
to two-point stencils 27 : 


cently approximated spatial derivatives to mixed space-time deriva- 
tives in a manner reminiscent of the Lax-Wendroff approach 36 by 
differentiating the following equations in space and time: 


duj 

~dT 


(m du > . 3 p\ d P ( M dp ■ a “A 
y ‘ ~dxi + 9xy J ’ ¥ = "( Mi a^ + 9^J 


(14) 


The mixed space-time derivatives are then used in modified series 
expansions with local coordinates about the center of the interpola- 
tion stencil at the current time level. For two-dimensional problems, 
the following series are used, with O = 2s 4- 1: 

o o 2(0) 

p(x, y, t) = cr ^yh' 

I 0 j 0*0 


O O 2(0) 

u(x, y, t) = J2 C “j,k x ‘y’ ,k 


i-o y - o * = o 


O 2(0) 


v(x, y , l ) = J2T,T. c ^ x ‘y 1,k 


i=0 j - o * - o 


(15) 


where 




/ j \ 3' --J 4 k f(x = 0, y = 0, t — 0) 
y i \j \k\ J dx 1 3 y } dt k 


d^f(x = 0) 
3X 1 * 1 


2S -r 1 

i - dr 


The variables p, w, and v and their spatial derivatives are then 
(11) advanced to the next time step by evaluating Eq. (15) and its spatial 
derivatives with x — 0, y = 0, and / = At as follows: 


where the function Z, ^ is independent of space and time and can 
be computed simply as follows: 

DO dx = 0, 2s + 1 


DO i — dx , s 


Z,,d* = 


(i — dx)! 


(-*o) 


0 ‘ -dr) 


END DO 


2 ( 0 ) 


^V<O.0.AO_ E(ali|)C , uai , 
* = 0 
2 ( 0 ) 


dx a y b 

3“ +(, d(0, 0, At) 


k = 0 
2(0) 


3 x a y b 


= 53(fl!6!)C; M At 1 


(16) 


i = 0 


END DO 

DOdx = 0, 2s + 1 
DO i = s + 1 , 2s + 1 


( 12 ) 


7 -i.Ax = T f r: — -TT , (-*o) (, ' t ' 1 ~ f) (-*i) 
^\(dr-r)!r! 


(i — s - 1 - dx + r) 


By unrolling the Cauchy-Kowalewski recursion in Eq. (14), we 
can quickly express the mixed space-time derivatives C% bk , C“ b k , 
and C v ab k in Eq. (16) in terms of the space derivatives C^ b 0 , C“ b 0 , 
and C* b 0 , which were approximated at the stencil center in the last 
section. The following loop efficiently unrolls the recursion 29 : 

DO* = 1,20 
DO 6=0, O 


dx-l-r r - 1 \ 

* D [* - (s + 1) - e] * ]> + 1 - i] ] 

e=0 k = 0 / 

END DO 

END DO (13) 

With these developments it is now possible to approximate effi- 
ciently all solution variables and their derivatives at the center of a 
two-point stencil to any level of accuracy with only a page of code. 
Because these explicit forms are completely algebraic, it is possi- 
ble to adapt dynamically the accuracy of the approximation to the 
solution evolution. 


DOa = 0,0 

a = a + 1, b = b+ 1, *=*-1 

C P a . b . k = {-{(a)*{M x *Cl fci + C“ fri )) 

= (-(<*> * M y * c« u ) - ( a> * (c; fc . + m x * c“ bi ))/k 
c:. b . k = (-((b) * (c; u + M y * c; u )) - (a) * M x * c; bi )/k 
END DO 


B. Time Evolution 

To utilize fully the arbitrarily high-order accuracy in space that is 
possible on two-point stencils, it is necessary to achieve similar accu- 
racy in time. This can be accomplished with the MESA method. This 
method uses the governing equations [Eqs. (1)] to convert the re- 


END DO 

END DO (17) 

This simple form applies to uniform mean flow, but for general 
flows Cauchy-Kowalewski recursion becomes complicated. 26 For 
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example, using Eq. (14) to calculate the mixed space-time derivative 
C“ 0 j , when the mean flow varies in space, results in the expression 


d 2 u / 3 Mi 3 u d 2 u 3 2 p\ 

3t3x y dx 3x, 1 3x ( 2 + dx 2 ) 


(18) 


which requires knowing the derivatives of the mean flow, dMi/dx, 
and the higher-order derivatives of Eq. (18) require higher-order 
derivatives of the mean flow as well. Getting this information may 
require computing the mean flow using a very high-order MESA 
method as well, and the equations will grow exponentially unless a 
simple form analogous to Eq. (1 1 ) is found. Despite these complica- 
tions, Cauchy-Kowalewski recursion, when automated with com- 
puter algebra tools, 29 can produce numerical schemes with higher 
order than is currently possible with Runge-Kutta methods. 


V. Including Wall Boundaries 

The procedure for advancing the interior grid points is not appli- 
cable to fill grid points near boundaries because part of their stencil 
is on or within a solid object. We obtain data for these points at each 
time step with a Hermite interpolant that is consistent with the wall 
boundary conditions and that uses the nearby solution at interior 
grid points. This Hermite interpolant is equivalent to the interpolant 
used for the interior grid points in Eq. (4) but is written as 


/(*. 


] s 

(*»(//,,. dv(y)) 

i,j 0 dx.d> ~ 0 


3* 4 dy /(*,>>) 

dx^dydy 


(19) 


where the H Xidx (x) and H yj dy (y) terms are 2s + 1 -order poly- 
nomials in x and y, respectively. Because the data at the interior 
grid points are known, this Hermite form of the interpolant has the 
advantage of reducing the number of unknown coefficients, C (- , 
in Eq. (4), to only the unknown data which are required at the fill 
points, 3 d * + d> / ( x t , y ; )/3x djr 3y d - v . 

Each shaded box in Fig. 5 represents an interpolation region used 
for interpolating the values at the fill points in each box. The bound- 
ary conditions are imposed upon the shape function for each box 
at the locations on the surface intersected by the arrows. The num- 
bers inside each box show the sequence in which the fill points are 
interpolated. 

A. Choice of Wall Boundary Conditions 

Each fill point on the grid provides (s + l) 2 unknowns in Eq. (19), 
so that an equivalent number of constraint equations must be ob- 
tained. An infinite number of constraints can be obtained by ex- 
ploiting that the physical boundary conditions apply uniformly in 
time everywhere in the boundary. 4 For Eqs. (1 ), the inviscid bound- 
ary condition is that the normal component of velocity is zero at a 
surface: 


Vrj = 0 


( 20 ) 


Because vorticity is convected with the mean flow, we also assume 
in this work that there is no vorticity at a surface: 


3v 

dx 



( 21 ) 


B B B 

B B B 

B B B 

B B B 

B B B 

B B B 

B B B 

B B B 

-B 4 * — — & 

-2 

B B B 

B B B 

B B B 

B B B 

B B B 

B B B 

B B B 

B B B 



B B B 
B B B 
B B B 


Fig. 5 Interpolating fill point data. 
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For flat walls, the boundary conditions for this case are 


3^- ] + T p 
3f} 2n + l 3x T 


d 2A + T V fl 
Sff^d i T 


= 0 , 


VnJ:nJe (0, 1,2,...) 


3 2« + i + r V{ 
drj 2n * l 3i T 


( 22 ) 


where V is the magnitude of the perturbation velocity vector with 
velocity components u and v in the Cartesian x- and y-axis direc- 
tions, p is the scalar pressure, fj is the unit vector normal to the wall 
surface direction, f is the unit vector tangent to the wall boundary 
direction, Vf is the velocity tangential to the wall, and V ^ is the 
velocity normal to the wall. 

Selecting which set of boundary conditions to use and where to 
apply them on the boundary can be simplified by considering the 
general form of the normal and tangential derivatives of f(x,y) 
(Ref. 37): 


a N + T f(x,y) 

dfj»dr T 






(23) 


If each operator is expressed in its binomial series form and the 
order of summation rearranged, then by substituting Eq. (19) into 
Eq. (23) the following form is derived: 


a N + T f(x, y ) ^ v' r v .a^/fa.yj) 

SfjNfT JL, 1 dz ’ d >' •*'•>'; dx^dyty 

1 ij — 0 dx.dv = 0 7 


(24) 


where Y^ x ^ y>Xit y i is defined by 

T N 

Ydx,dy.x it yj = K( a > b) 


d N + T (H Xitdx (x)H yjA y (y)) 

3x N +T o-bfya+b 


(25) 


K( °' )_ ((JV-Z>)!6!(7'-a)!a! / )' ,J[ 


b r\l a+b (-\) T “ 


(26) 


Because maintaining numerical stability for hyperbolic problems 
normally requires methods that include all information from within 
the domain of dependence defined by the characteristic surfaces, 38 
we want to select N = 2n + 1 and T so that 

3x N ~ T -a - bftya + b 


in Eq. (25) never becomes zero because this would exclude 
some of the influence of neighboring grid points. For example, if 
Y&xAy'Xi.y, =0, then the grid point data y J )/3x dx 3y dy 

would be excluded from the boundary condition in Eq. (24). 
Similarly, if a part of the sum forming Y^ x ^ y ^ yj is missing, 
then some influence of the corresponding grid point is absent as 
well. Maintaining the full influence is accomplished by choos- 
ing (n :n =0, 1,2 5 ) and [T : T — 0, 1, 2, . . . , 2(5 -n)]forthe 

first and third boundary conditions in Eq. (22). Those conditions 
will produce linearly consistent systems that can be solved to de- 
termine the interpolant or, equivalently, the data required at the 
fill point. However, in the special case in which rj x = 0 or r) y = 0, 
we need to use the different conditions (n : n = 0, 1 , 2, . . . , j) and 
(T : T = 0, 1, 2, . . . , s) to ensure nonsingular matrix systems. This 
special case does not result in a loss of information because the 
affected terms in Eq. (25) are supposed to be zero. For example, a 
third-order MESA method would use the following pressure bound- 
ary conditions at each location indicated by the arrows in Fig. 5: 



d 2 P 

df)dx 


= 0, 


dP 

3fjdi 2 


= 0, 



(27) 


The same method applied to the special case shown in Fig. 6 would 
use the conditions 



d 2 p 

3rjdz 


= 0 , 




(28) 


B. Symbolic Solution 

The shape function for each fill point can now be solved using 
the boundary conditions and the known interior data in each shaded 
region of Fig. 5. The boundary condition Eqs. (22) are applied to 
the shape function Eq. (24) to solve for the unknown data at the fill 
points. This linear system of (s + l) 2 equations per fill point may 
be written as 


Mf = Nd (29) 

where M and N are matrices whose elements are given by ldx.dy,*,.^ 
in Eq. (25), where / is a vector of the data needed for the fill points 
and where d is a vector of known data from the interior points. 

This system can be solved using any linear solver at each time step, 
but for high-order methods this becomes expensive and the matrices 
become ill conditioned. A better approach is to solve symbolically 
this system once to form the algebraic solution 

/ = M ] Nd (30) 

where M and N are numerical matrices dependent on the geometry 
but the vectors are symbolic. This results in an algebraic expres- 
sion for each fill point as a linear combination of the interior grid 
point data. It is then a simple matter to evaluate the fill points by 
updating the vector d at each time step and evaluating their linear 
combinations. 

This approach worked well for up to 1 lth-order accuracy in two 
dimensions for the case (shown in Fig. 6) of no mean convection in 
an unrotated box that is aligned with the coordinate axes. For higher- 
order algorithms, the matrix M becomes ill conditioned, and finding 
its numerical inverse becomes difficult, inefficient, and unstable. 
Resorting to Gaussian elimination would be more stable for high- 
order cases, but the cost at each time step is prohibitive. 

VI. Stencil Constraint Tree 

Instead of using a poorly conditioned numeric matrix M, its sym- 
bolic form can be inverted. There are actually only a relatively few 
symbolic matrices that need to be inverted regardless of the ge- 
ometry. These few cases can be found by using a new tree data 
structure, referred to as a stencil constraint tree (SCT), that repre- 
sents all possible stencil configurations. For example, in Fig. 7, all 
possible 3x3 stencil configurations are shown in which 1 ) the top, 
center grid point is a boundary point (indicated by B); 2) the center, 
center grid point is a fill point (open circles); and 3) the bottom, 
center grid point is an interior point (filled-in circles). The SCT for 
the simple case of a 2 x 2 stencil that has an interior grid point in its 
top-right location is shown in Fig. 8. It is constructed by propagat- 
ing symbolic constraints in a manner analogous to Waltz’s symbolic 
constraint propagation algorithm (see Ref. 39). 

Once the tree is constructed, a particular stencil configuration is 
found by traversing the SCT from top to bottom along any path and 
converting the branch labels using Table 1 and Fig. 9. For example, 
the branch number 10 corresponds to position 4 and label 1 in Table 1 
and represents an interior grid point that is located at the top-right 
comer of the stencil. A similar procedure is applicable to larger 
stencils as well. Notice that identifying a particular n -point stencil 
using an SCT can be done in n 2 or n 3 steps for two- and three- 
dimensional stencils, respectively, instead of the 3" n 2 and 3" n 3 
comparisons typically required using a brute-force comparison of 
all stencil configurations. 


Table 1 Branch number legend 


Position 

Label 

Branch 

1 

1-3 

1-3 

2 

1-3 

4-6 

3 

1-3 

7-9 

4 

1-3 

10-12 
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Fig. 7 All possible S7 stencil configurations with fill at center: 16 cases. 


A. Constructing the Tree 

The most significant advantage of the SCT is it can reduce the 
number of symbolic matrices from Eq. (29) that need to be consid- 
ered. This is accomplished by propagating natural constraints during 
the construction of the tree. 29 Natural constraints are simply a list of 
rules that are a natural consequence of the definitions of interior, fill, 
and boundary grid points on a Cartesian grid. For example, some 
natural constraints that limit the set of stencil configurations are 
1) an interior grid point is never adjacent to a boundary grid point 
and its converse and 2) for two- and three-point stencils, a fill grid 
point will always be adjacent to at least one boundary point. If the 
natural constraints are not enforced, then the list of stencil configu- 
rations in Fig. 7 grows to 3 6 = 729 because only the middle column 
is defined. The advantages of enforcing natural constraints increase 
with both spatial dimension and stencil size. 



Fig. 8 Two- point stencil constraint tree. 
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Fig. 9 Stencil grid positions. 
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Fig. 10 Sub 3x3 stencil sym- 
metry. 


S8 


S7 


B. Small Stencils Guarantee CFL Stability 

These techniques can be applied to both of the 5 x 5 stencils 
shown in Fig. 10 to construct their SCT. This SCT will represent all 
possible 5x5 stencils that can occur with a fill point at its center. 
We find that there are only 65 unique rotationally symmetric cases 
of 3 x 3 stencils centered about the fill point. 29 All other possible 
stencil configurations are simply rotations of this core group. Note 
that we assume there will be at least one interior point next to the 
fill point, otherwise that fill point is not needed as shown in the 
comers of Fig. 5. As shown in Fig. 10, there will always be a 3 x 3 
stencil contained within the 5 x 5 stencil that does not contain a 
boundary point and yet always contains a fill point on its edge or 
comer. This guarantees the existence of a 3 x 3 stencil for each fill 
point in any geometry that is not intersected by the surface geometry. 
Therefore, it will always be possible to expand the stencil outward 
from each fill point to the surface where the boundary conditions are 
imposed. This in turn guarantees that the CFL stability criterion will 
be satisfied because the domain of dependence for the fill interpolant 
is increased from a standard stencil size. 

This CFL guarantee is a property of two-point and three-point 
stencils in two or three dimensions. Four- point stencils or larger will 
require smaller time steps near the boundary and make it difficult to 
choose locations on the surface that guarantee a linearly consistent 
system in Eq. (29). 


= Possibly Degenerate 

Fig. 11 S8 symmetrical mapping: 18 cases. 
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Fig. 12 Use regions 2 and 3 instead of region 1 . 


A. Efficient Retrieval of Fill Point Solutions 

Once all of the stencil configurations are known and the direction 
in which to map each fill point to the surface is found, then the 
symbolic form is known for every matrix M and N in Eq. (29). 
Each of these systems may then be symbolically solved, and their 
subsequent solutions may be stored in the leaves of the SCT for fast 
retrieval later. 

The leaves of the tree will contain algebraic solutions for every 
fill point in any geometrical configuration. They can be used and 
reused many times without the need for solving Eq. (29) again. For 
example, if the walls are moving then the Cartesian grid points will 
need to be relabeled, but the same SCT, mapping, and symbolic 
solutions can be used at each time step to interpolate efficiently the 
data needed at each fill point. 


VII. Mapping 

All possible 5x5 stencil configurations containing a fill point are 
known after the construction of the SCT for the stencils in Fig. 10. 
This information is used to select locations on the boundary at which 
the boundary conditions are enforced. A surface will always occur 
between a fill point and a boundary point, so that even before we 
know exactly where the surface is, we know in which direction to 
go to find the surface. 

We would like to use the points on the surface that are closest 
to the interpolation stencil while ensuring that never more than N 
locations are collinear with the coordinate axes for each A -point 
interpolation stencil. These conditions maximize the accuracy of the 
interpolant and ensure that the matrix M in Eq. (29) is nonsingular. 
One such mapping is shown in Fig. 1 1 for the upper triangular part of 
the outlined S8 3 x 3 stencil in Fig. 10. The degenerate cases occur 
because it is not guaranteed that a nearby surface will be found using 
the mapping shown. These cases are usually not an issue in practice 
and can be avoided by choosing alternate interpolation regions as 
in Fig. 12 or by constructing the SCT for a larger stencil around the 
fill point in Fig. 10 and then choose a mapping with this additional 
information. 

Implementing higher accuracy boundary conditions with this ap- 
proach requires increasing the depth of Hermitian data at each grid 
point and increasing the number of boundary conditions used from 
Eq. (22), but the stencil size used for interpolation is never changed. 
Therefore, we never need to remap the fill points to the boundary 
and the CFL stability is not affected. 


VIII. Numerical Results 

All of the following results solve the linearized Euler equations 
[Eqs. (1)] in uniform mean flow using the approaches already de- 
scribed. This approach to developing very -high order methods was 
applied to both wave propagation and acoustic scattering problems 
in various test geometries. 

A. Acoustic Propagation 

For mean convection velocity M in an open biperiodic (triperi- 
odic) domain (for d = 2 or 3 dimensions), the linearized Euler equa- 
tions have the following analytical solution: 

d d 

/>(*, 0 = c, n V HfCr. ') = -S,.iC x , LI S x> (31) 

i — 1 j t i 

7 = 1 


where 


C Xl — cos(W t 7iXj), S Xi = sin(Wi-n’ii) 


C t =cos(jr/||W||), 


S,.i = 


W, sin(7r? || W||) 

\\w\\ 


IIWII = 



x, = (Xj - M,t) 
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Fig. 13 Resolution, efficiency, and memory usage of higher-order 
methods. 


The properties of very high-order methods were tested on 
this problem with mean velocity Af = (l, 1) and wave number 
W = (1. 1), and the results were used to validate the automatically 
generated codes. The methods were compared for efficiency and 
resolution; the results for up to 29th-order accuracy (using 64-bit 
precision) are shown in Fig. 13. Notice that higher-order meth- 
ods perform better in each part of Fig. 13 until machine preci- 
sion is reached. By increasing the wave number in Eq. (31) to 
W — (Wi > 2, W 2 > 2), we also see that very high-order methods 
offer subgrid scale resolution as shown to the left of the ordinate 
axis at the top of Fig. 13. Similar results were found while solving 
three-dimensional cases as well. 29 

The resolution of the methods improves while using 128-bit pre- 
cision as the accuracy increases until about the 21st order for signals 
with eight grid points per wavelength: This limit varies with signal 
frequency. If the order is increased above this limit, then the resolu- 
tion of each signal degrades from roundoff errors that begin to dom- 
inate around 57th-order accuracy for signals with a wavelength of 
eight grid points. However, by selecting wave number W = (16, 16) 
in Eq. (3 1 ), and by using 40th-60th order methods on a two-point 
stencil, it is possible to resolve 4 wavelengths per grid point. Note 
that this is not a violation of the Nyquist criterion 40 because each 
grid point contains multiple data values for Hermite interpolation. 
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Fig. 14 Resolution of rotated boxes and circles. 
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Fig. 15 Unrotated box grid resolution studies. 


For extremely high-order methods, the convergence slope at the 
top of Fig. 13 is so steep that essentially only a specific frequency 
is correctly resolved with finite precision floating point hardware 
unless a method for detecting roundoff error is available. Fortu- 
nately, roundoff error may be detected by checking the columns of 
the Hermite divided-difference tableau in Fig. 4, which is used for 
the interpolation. The sum of column j should equal the difference 
of the first and last terms in column j - 1 (Refs. 27 and 41). If 
roundoff error begins to occur, simply stop constructing the tableau 
at that point and use the resulting lower-order method. This saves 
computational effort and enables the solution of a field containing 
widely varying wavelengths. 


B. Acoustic Scattering 

With the exceptional qualities of very high-order methods estab- 
lished for wave propagation, the simulation of acoustic scattering 
from within a rotated box and a circle was completed (Fig. 14). An 
analytical solution to wave scattering within the rotated box with no 
mean convection is 

p(x, y, /) = - cos(v / 27rO cos(7rx) cos(7ry) (32) 

u( jc, y, r) = -[cos(7ry) sin(V / 27r/) sin(7rx)/V2] (33) 

v(x, y, 0 = — [cos(7rx) sin(V^7r/) sin(7ry)/V21 (34) 


where the coordinate system is aligned with the box. An analytical 
solution to wave scattering within a circle with no mean convection is 


5J 0 (kr) cos(A-f) 


(35) 


where A. = 3.83171, Jo is the Bessel function of the first kind of 
order 0, and polar coordinates are used. 

For the case of an unrotated box, we found stable solutions 
with 2nd-, 3rd-, 5th-, 7th-, 9th-, and llth-order accuracy using 
boundary conditions from Eq. (22) with (n : n = 0, 1 , 2, . . . , s) and 

( T :T — 0, 1, 2 j) as shown in Fig. 15. Using the boundary 

conditions in Eq. (22) with n, T — 0, we found stable second-order 
solutions for all box and circle cases, regardless of the rotation an- 
gle. The first five plots in Fig. 14 correspond to box rotation angles 
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0, n/ 16, 7T/8 , tt/ 4, and 7r/3, respectively. The last curve in Fig. 14 
shows a second-order method without wall boundaries. 

C. Numerical Stability 

For second-order boundary treatments, ordering interpolation re- 
gions by the number of fill points they contain from fewest to most 
and then interpolating the fill points in that order resulted in sta- 
ble solutions for irregular geometries. There are, however, multiple 
stable choices of shaded regions and interpolation sequences, and 
rather than relying on numerical experiments, an important next step 
is to use the symbolic form in Eq. (29), of which only a relatively 
few cases are possible, to perform a complete stability analysis. 
With this better understanding, perhaps higher-order conditions can 
be made reliably stable in irregular geometries by explicitly adding 
artificial dissipation or by rearranging the interpolation regions. 

IX. Conclusions 

Small stencils containing Hermitian data were found to possess 
key properties necessary for developing exceptionally efficient, ex- 
plicit algorithms of arbitrary formal accuracy in space and time 
for scattering problems with arbitrary surface orientations. This is 
accomplished by approximating spatial derivatives using a special 
form of two-point Hermitian divided -difference spatial interpolation 
and an unrolled Cauchy-Kowalewski recursion procedure for time 
advancement. A stencil constraint tree is used to find all possible 
stencil configurations and their respectively well-posed boundary 
conditions so that grid points near arbitrarily oriented surfaces may 
be interpolated. Computer algebra is then used to find the symbolic 
shape function for each stencil configuration, and the stencil con- 
straint tree is used to identify quickly the correct interpolant for each 
near boundary grid point. 

In spite of the development of simple procedures for the time evo- 
lution at all grid points, the FORTRAN application code of a very 
high-order method is complicated, and so the task of programming 
was automated using computer algebra procedures. These proce- 
dures have successfully produced 1 ) stable, parallel 57th-order wave 
propagation methods in two and three spatial dimensions with peri- 
odic boundary conditions, 2) stable 2nd-order methods for scatter- 
ing problems with generalized surface boundary conditions in two 
spatial dimensions, and 3) stable 1 lth-order methods for scattering 
problems with surfaces aligned to the grid in two spatial dimensions. 

The resolution and efficiency of these methods improve with or- 
der until the accuracy exceeds the limits of machine precision. At 
these very high orders of design accuracy, it is necessary to control 
roundoff error, and a method for doing this is suggested. 

The procedures discussed here provide a systematic method for 
quickly changing the accuracy of both the interior propagation and 
wall boundary condition implementations to suit flow conditions. It 
was shown that very high-order methods (>15) provide an oppor- 
tunity for subgrid-scale resolution. 

These results demonstrate the potential value of very high accu- 
racy in time as well as in space for aeroacoustic calculations. How- 
ever, a detailed stability analysis for these procedures remains to be 
done for the generalized high-order surface conditions. A method 
for unrolling nonlinear Cauchy-Kowalewski recursion needs to be 
developed to achieve arbitrarily high accuracy in time for nonlinear 
problems. Finally, arbitrary order accuracy nonreflecting boundary 
conditions need to be developed to take full advantage of the arbi- 
trary interior domain accuracy. 

Nonetheless, it is necessary to use Hermitian data for achiev- 
ing very high-order solutions because Lagrangian stencils simply 
become too large. Because since small stencils allow for a very 
simple interpolation procedure of arbitrary accuracy while main- 
taining CFL stability near boundaries, future work will be devoted 
to exploiting them in substantial applications. 
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